Constraint preserving boundary conditions for the Z4c formulation of general 

relativity 
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We discuss high order absorbing constraint preserving boundary conditions for the Z4c formu- 
lation of general relativity coupled to the moving puncture family of gauges. We are primarily 
concerned with the constraint preservation and absorption properties of these conditions. In the 
frozen coefficient approximation, with an appropriate first order pseudo-differential reduction, we 
show that the constraint subsystem is boundary stable on a four dimensional compact manifold. 
We analyze the remainder of the initial boundary value problem for a spherical reduction of the Z4c 
formulation with a particular choice of the puncture gauge. Numerical evidence for the efficacy of 
the conditions is presented in spherical symmetry. 
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I. INTRODUCTION 

Numerical simulations of general relativity typically 
introduce an artificial time-like outer boundary. This 
boundary requires conditions which ought to render the 
initial boundary value problem (IBVP) well-posed. Well- 
posedness is the requirement that a unique solution of the 
IBVP exists and depends continuously upon given initial 
and boundary data. 

The most important approaches to demonstrate well- 
posedness of the IBVP are the energy and Laplace- 
Fourier transform methods [H-Q • The energy method is 
straightforward. In this approach one constructs a suit- 
able norm for the solutions of the dynamical system. Us- 
ing the equations of motion one can estimate the growth 
of this norm in time. However, this technique in general 
cannot be used if the system is not symmetric hyper- 
bolic or if the boundary conditions are not maximally 
dissipative. Recently, Kreiss et al. introduced in Q a 
non-standard energy norm to prove that the IBVP for the 
second order systems of wave equations with Sommerfeld- 
like boundary conditions is well-posed. The key idea 
in [f| is to choose a particular time-like direction in a 
way that the boundary conditions are maximally dissi- 
pative ones. A different method is based on the frozen 
coefficient approximation. In this approach one freezes 
the coefficients of the equations of motion and the bound- 
ary operators. The IBVP is thus simplified to a linear, 
constant coefficient problem which can be solved using a 
Laplace-Fourier transformation. Sufficient conditions for 
the well-posedness of the frozen coefficient problem were 
developed by Kreiss in Q if the system is strictly hy- 
perbolic. Using that theory, a smooth symmetrizer can 
be constructed with which well-posedness can be shown 
using an energy estimate in the frequency domain. Agra- 
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novich [7| extended that theory to the case in which the 
system is strongly hyperbolic and the eigenvalues have 
constant multiplicity. It is expected that, by using the 
theory of pseudo-differential operators Q, one can show 
well-posedness of the general problem. In what follows 
we use the Laplace-Fourier approach to prove the well- 
posedness of the IBVP for the constraint subsystem of 
Z4 with high order constraint preserving boundary con- 
ditions. The order of a boundary condition refers to the 
highest derivative of the metric or of the gauge variables 
contained therein. 

Once the continuum boundary conditions (BCs) are 
understood, one needs a strategy for their implementa- 
tion in a numerical code. The numerical implementation 
is required to be stable and to converge to the contin- 
uum solution at a certain rate. If the system is sym- 
metric hyperbolic one can use, for instance, difference 
operators which satisfy summation by parts schemes and 
penalty techniques to transfer information through the 
outer boundary condition [9l4T2l|. This allows the deriva- 
tion of semi-discrete energy estimates which can guaran- 
tee the stability of the numerical implementation. Nev- 
ertheless, in the general cases, even for a linear system, 
demonstrating that a numerical approach to BCs will re- 
sult in a stable scheme is difficult. In the absence of a 
proof of numerical stability one must rely on calculations 
for similar toy problems and on thorough numerical tests 
with simple data, in which problems can be identified 
locally at the boundary. Unfortunately naive discrete 
approaches are often numerically unstable. One issue is 
that a code usually requires more conditions than are 
given at the continuum level. 

The two most popular choices of formulation of general 
relativity (GR, hereafter) in use in numerical relativity 
today are the generalized harmonic gauge (GHG) and 
the BSSN formulations EMI. Significant progress has 
been made in the construction of both continuum and 
discrete BCs for the GHG formulation, see e.g. fl7H24j ] 
and references therein. For GHG the task is made rela- 
tively easy because the system has a very simple wave- 
equation structure in the principal part and furthermore 
because the constraints may be expressed as time deriva- 
tives of metric fields. The BSSN formulation is used to 
evolve both vacuum and matter space-times by a num- 
ber of numerical relativity groups, see e.g. |25l - [32j and 
references therein. This The system is taken in combina- 
tion with the so-called moving puncture-gauge [133 • BCs 
for the BSSN formulation have received relatively little 
attention, although recently, Nunez and Sarbach have 
proposed [331 ] constraint preserving boundary conditions 
(CPBCs) for this system. Recasting the dynamical sys- 
tem into a first order symmetric hyperbolic system, they 
are able to prove that the corresponding IBVP is well- 
posed through a standard energy method, at least in the 
linearized case. These boundary conditions have not yet 
been implemented numerically. Currently Sommerfeld 
BCs are the most common in use in applications, despite 
the fact that they are certainly not constraint preserv- 



ing and it is not known whether or not they result in 
a well-posed IBVP. The problem is that the character- 
istic structure of puncture-gauge BSSN is more compli- 
cated than that of the GHG formulation, which makes 
the analysis difficult. Despite the fact that with Sommer- 
feld conditions the constraints do not properly converge, 
in applications they are robust and are currently not the 
dominant source of error in numerical simulations. 

Another version of GR is the Z4 formulation [U, [35[ . 
When coupled to the generalized harmonic gauge Z4 is 
formally equivalent to GHG [36| . Additionally it is possi- 
ble to recover the BSSN formulation from Z4 by freezing 
one of the constraint variables. In this sense Z4 may 
be thought of as a generalization of both BSSN and 
GHG. Z4 has the advantage over GHG that it main- 
tains sufficient gauge freedom that it may be coupled to 
the puncture-gauge, potentially allowing the evolution 
of puncture initial data as is standard with BSSN. To 
that end a conformal decomposition of the Z4 formula- 
tion (hereafter Z4c) was recently presented [37j . Unlike 
BSSN, the Z4 formulation has a constraint subsystem in 
which every constraint propagates with the speed of light, 
which may be useful in avoiding constraint violation in 
numerical applications. It was shown that, at least in 
the context of spherical symmetry, numerical simulations 
performed with puncture-gauge Z4c have smaller errors 
than those performed with BSSN [37j . However it was 
also found that Z4c is rather more sensitive than BSSN 
with regards to boundary conditions. BCs compatible 
with the constraints for a symmetric hyperbolic first or- 
der reduction of Z4 were specified and tested numerically 
in [38|, . Those conditions are of the maximally dis- 
sipative type and, therefore, the well-posedness of the 
resulting IBVP was shown by using a standard energy 
estimation. Nevertheless, Bona et al. used in [38l |39| 
harmonic slicing and normal coordinates to rewrite Z4 
as a symmetric hyperbolic formulation. Therefore, it is 
not clear if their results can be easily extended to the 
general case. 

In this work we therefore specify BCs in combination 
with puncture-gauge Z4c and we show that the resulting 
IBVP is well-posed at least in spherical scenarios. How- 
ever, since we are interested in specifying CPBCs which 
can be used in 3D evolutions, the well-posedness of IBVP 
for the constraint subsystem is established in the general 
case. In addition, we study the effectiveness of these con- 
ditions by performing numerical evolutions in spherical 
symmetry. 

We begin in section|H]with a summary of the Z4 formu- 
lation, and identify the BCs we would like to consider in 
our analysis. We present the analytic setup for our well- 
posedness analysis in section Hill Section ITVl contains our 
analytic results on BCs for Z4. We present our numeri- 
cal results in section IVl We conclude in section IVTl The 
principal ideas of the Kreiss theory are summarized in 
appendix [X] and applied to the wave equation with high 
order BCs in appendix [Bj We describe the numerical im- 
plementation of the second-order CPBCs in Appendix [Cl 
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II. THE Z4C FORMULATION 

In this section we present the Z4c formulation and the 
general expressions for our BCs in detail. We also intro- 
duce notation for a 2 + 1 decomposition in space, which 
we will use in the calculations in the following sections. 



properties of the formulation. In what follows we will use 
the shorthand 

f d * = f^kl~. kj = 7 i 7 ii 7 « ^ _ ^_ lu ^ (13) 

In terms of the conformal variables the evolution equa- 
tions for the Z4c formulation become 



A. Evolution equations and constraints 

Following 37], in which the conformal Z4 formulation 
was presented, we replace the Einstein equations with the 
expanded set of equations 



2 a K r 



(1) 



DiDja 



Kij K 



ae = a 



a [Rij-2K ik K, 
+ 2 dyZfi] +Aira [ 7y (S-p)-2 S y ] 
+ CpKij , (2) 
1 



H + d k Z k 



p e, 



d t Zi = aM i + a@, l + P 3 Z,, 



(3) 
(4) 



where 9 and Zi are constraints. The ADM equations are 
recovered when the constraints and Zi vanish. The 
Hamiltonian and momentum constraints 



H = R- K lJ + K 2 - 16 up = , 
M 1 = Dj (K ij - j ij K) - 8 7T S 1 = , 

evolve according to 



dnH 



-2diM\ 
1 



(5) 
(6) 

(7) 
(8) 

(9) 



doM, ~ -- d t H + d 3 d,z t - di&Zj , 

in the principal part, where we have used 

flb = - (dt - ^flfc) . 
a 

Since we are concerned in this work only with the behav- 
ior of the BCs on the constraints, we have discarded the 
constraint damping scheme of |36l |. 



B. Conformal decomposition 

In our numerical applications we evolve the Z4c system 
in the conformal variables Xi %j , K , Aij and T l , defined 
by 

7« = l~ k 1*3 . K = YSKu -29, (10) 

X = 7"' , Aj = 7"^ {Kij - \lijK) , (11) 
and finally 

r i = 2f j Z j +f^% k ,i. (12) 

The choice of conformal variables allows us to evolve 
puncture initial data whilst altering the underlying PDE 



dtX = %XHK + 26) - Dip] , 



(14) 



Octo = -2aAij + p% ik + 2j {ilk p« {j) - -%P*,k , (15) 

d t K = -D l D t a + a[A tJ A ij + -(K + 26) 2 ] 

+ 4na[S + Padm ] + (3 l K, t , (16) 
the trace-free extrinsic curvature evolves with 
d t Aij = xl-D.Dja + a(Rij - ^irS^f 
+ a[{K + 2e)Aij-2A k iA k j] 

+ (3 k A^ k + 2A [l]k p\ b) - ^A tJ /3 k . k , (17) 



and finally we have 



3 r, 



d t V = -2A^a, 3 + 2a[T) k A 3k - - A 13 ln( X ) 



r(2K + e),j-8^S j }+j 3k p% 



Z&3 

3 
1 

3 



31 ' / rjk 

2 - 



+ ff ij P k kj + PKj - f d% + §IY/% . (18) 

The <d variable evolves according to Eqn. §5§ with the 
appropriate substitutions Eqns. (|22II25[) . In the dtAij 
equations we write 

R i j = RX is + Rlj t (19) 
^DaDa - ^ 3 D l xD lX , (20) 



Rij = --^l lm iij,lm +7fc(i|f l | )J -) + fd f(y)fc + 

7 " n (2f? (j f J - )fcro +f? m f Wi ) . (21) 
The complete set of constraints are given by 

6 , 2Z % = jijt 3 - yi%j, k , (22) 

H = R- A»Aij + \{K + 29) 2 - 16tt Padm , (23) 
IVP = djA ij + r*j k A ik - p ij dj(K + 26) 

_|i*i(logx),i, (24) 

D = ln(det 7) = , T = f 3 A VJ = . (25) 

In our numerical evolutions the algebraic constraints D 
and T are imposed continuously during the numerical 
calculations. It seems to improve the stability of the 
simulations significantly [40l |. 
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C. Puncture gauge conditions 

The most popular gauge choice in the numerical evo- 
lution of dynamical spacetimes is the puncture gauge. In 
introducing scalar functions (fi^, jj,s,e a ,e x ) the general 
form of the gauge (without introducing the additional 
field B l ) is 

dtd = fru.i — (1l a 2 K , (26) 

+ e x f j d jX - (27) 

Note that in this condition we have included a new term 
proportional to the spatial derivative of X- As we wu l 
show in Sec. IIVI the inclusion of that term allows us 
to prove the well-posedness of the IBVP for Z4 in the 
frozen coefficient approximation. We refer to the choice 
of shift (^s,e a ,e x ) = (1,1,1/2) as the asymptotically 
harmonic shift condition because in preferred coordinates 
in asymptotically flat spacetimes near infinity the condi- 
tion asymptotes to the harmonic shift. 

In evolutions of equal mass black holes the gauge 
damping parameter r\ is usually taken as 2/M, where 
M is the ADM mass of the spacetime. Recently it has 
been shown that a spatially varying r\ parameter ma y b e 
helpful in the evolution of unequal mass binaries [4l],|42j ■ 



D. Fully second order form 



The principal parts of the constraint subsystem are just 
wave equations 



□9 ~ , UZi ~ , (33) 
DH ~ , DM, ~ . (34) 



Following the approach of [43j we will analyze the system 
starting in fully second order form. The equations of 
motion can be 2 + 1 decomposed against the spatial unit 
vector s l . We define the projection operator q l j — S l j — 
s l Sj. The equations of motion split into scalar, vector 
and tensor parts. The decomposed variables are written 

7 SS = sV7y , j gg = q l: ><y i: j , (35) 
J sA = sV Alij , (36) 

Tab = ~ \ub^ lij , (37) 

p s = s l p i , p A = q l A Pi • (38) 

The metric and shift are reconstructed from the decom- 
posed quantities by 



The principal part of the Z4c formulation in fully sec- 
ond order form is given by 



{d 2 a - fiLdid 1 ) a ~ . 



(28) 



7 3 Ms 
a 2 n L 



e a d d l a 



1 / 7 3 Ms 
3a I 2 



l^dod^jk, (29) 



(a 2 - d t d l ) 7 „ * V ( i - 4^) 

6 V 7 3 MS/ 

didjot (30) 
lk(td 1 )d p k . 



1 



7 3 Ms, 
a 2 e c 



" v 7 3 Ms 
or 



- 1 



One may view the constraints and Zi as being defined 
by the gauge choice (126H27p . 

2 6 = — d Q a - W j d olij + -d t f3 l , (31) 
aHL 2 a 



2Z, 



Ms 7 



(a jij d p j +i]Pi + e a a a :i 

-e x7 1/3 a iX )-(f d ). i . (32) 



(A R 1 j d \ TP A 

SiSjjss + q tJ j qqi (39) 
Pi = Sip s + qfp A . (40) 



Here and in what follows we use upper case Latin indices 
to denote projected quantities. 



E. Characteristic variables 



The standard parameters choice in the gauge condi- 
tions is fiL — 2/a, fis — 3/4 and e a = e x = 0. When Z4 
is coupled to the puncture gauge it is typically strongly 
hyperbolic (necessary and sufficient for well-posedness of 
the initial value problem, see e.g. |44j) except in a 
handful of special cases which for brevity we do not dis- 
cuss here. The fully second order characteristic variables 
with e a = e x = were presented in [32} ■ Here we present 
them with the additional parameters. In the scalar sector 
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the characteristic variables are 



U 



= d a ± y/JIZdsa , 



(41) 



U ±x = d A ± " A d s A 

7 ' Ms 

a 2 -7 1/3 /is 



a A 



— 5 ^ -p: r— dsfis ± — TJT, d fi. 

(A or — I 1 ' J /is) OL \ 7 ' /is 



2 (A a 2 +7 7 Mi Ms) 



(A 2 — /xl)(A 2 a 2 — 7 1 / 3 /iL /is) "Mi 
2(e a - l)aA 



a /iLAe a 
7 1/3 /is 



<9 s a 



(A 2 - AtL )(A 2 a 2 - 7 1 /3 Ms ) 

2 (A 2 7 1/3 /is + Q 2 e af j, L ) 
(A 2 — /^l)(A 2 a 2 — 7 1 / 3 /is) a/iL 



(<9 s a =F A 9 «) 



woQ ± — tt^ o s a 

7 ' Ms 



1/ = 907>J9 ± ds^qq , 

,+1 7 1/3 -Y 1/3 WS 

a /i_L 



7 1/3 MS 



2 a 



3a 



0.A. 



(42) 
(43) 



(44) 



with geometric speeds ±( y /jr[ 1 A, 1, 1) respectively, and 
where we have defined 



A = 7 SS + 7gg , 



A 



/ 2 (2 7 i/3 Ms -e x ) 
3 a 2 



(45) 
(46) 



Note that if A vanishes the system is only weakly hy- 
perbolic. In the vector sector the characteristic variables 
are 

vV3 



= d f3 A ± VMS7 9 s /3^ , (47) 



1 = 5 o7sA ± dslsA 1/3 



7 Ms 



-(a /3A±a s /3 A ), 



(48) 



with speeds ±( V //Is, 1). Finally in the tensor sector we 
have simply 

U±l = d 01 Z±d sl H, (49) 

with speeds ±1. 

Since the constraints (9, Zi, H, M,-) satisfy wave equa- 
tions in the principal part their characteristic variables 
are simply 

U±=d O±d s O, Uf = d Z s ± d s Z s , (50) 

U±=d Z A ±d s Z A , U±=d H±d s H, (51) 

= d M s ± d s M s , U% = d M A ± S M A , (52) 
each with speeds ±1. 



F. High order absorbing constraint preserving 
boundary conditions 

Following the notation of [23| we investigate the Z4c 
evolution equations on a manifold M = [0, T] x S. 



The three dimensional compact manifold S has smooth 
boundary 9E. The boundary of the full manifold T = 
[0, T] x <9£ is timelike and the three dimensional slices 
S f = {<} x S are spacelike. The boundary of a spatial 
slice is denoted St — {t} x dT,. 

We define a background metric (&, /3j, 7^ ), 



ds 2 = </ ab dir a da; 6 = 
- a 2 dt 2 +% j (dx i 



/3 i dt)(dx j +p j dt). (53) 



We assume the background 3-metric to be conformally 
flat for later convenience. It can be written as 



7y dx* dx j = ip 4 (dr 2 + r 2 dn 2 ) , 



(54) 



which defines the background isotropic radial coordinate 
r and the metric on the two-sphere is dSl 2 . We further- 
more define h a , the background future pointing unit nor- 
mal to the slices £(, and s a , the unit background nor- 
mal to the two-surface {t} x 9E as embedded in £ t . 
We are primarily concerned with absorbing conditions 
for the Z4c formulation as a PDE system. Construct- 
ing BCs explicitly related to the incoming gravitational 
radiation is left for future work. For simplicity our time- 
like and outgoing normal vectors (h a ,s a ) are therefore 
defined against the background metric. Therefore, 



9abn n 



1 



9abS a S b 



1 



g ab n a s b = 0. (55) 



To finish formulating the BCs we define the back- 
ground outgoing characteristic vectors 



(56) 
(57) 
(58) 
(59) 



where v s and vt are the characteristic speeds associated 
with equation (12T)1) in the scalar and vector sector, respec- 
tively. Since the constraints 0, Zi satisfy wave equations, 
constraint preserving, absorbing BCs in the linear regime 
around the background are given by 



(r 2 l a d a ) 6 = 0, (r 2 l a d a y 



Zi = 0, 



(60) 



where L > is an integer and = denotes equality in 
the boundary T. Note that the above conditions can be 
considered a generalization those recently proposed by 
Bona et.al. in [38|, |39( which correspond to L = (see 
also ©HI). 

We assume that both the physical and background 
metrics are sufficiently close to flat so that the full system 
has ten incoming characteristic variables at the bound- 
ary, which determines the number of BCs we may specify. 
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The boundary conditions (jST))) give four of the total. We 
take for the remainder 



,2 Juap, \ L + 1 


— h 




\ L+l 

2 k a d a j p s 


= h s , 


(62) 


\L + l 

2 J a da) Pa 


= h A , 


(63) 


1 Oa) Tab 


— h TF 

— n AB ' 


(64) 



where h a ,hi , h\^ B are given boundary data. Since the 
above conditions are not tailored to the characteristic 
structure of the puncture gauge Z4c system there is 
no guarantee that they will absorb all outgoing fields. 
We leave detailed examination of the absorption proper- 
ties, and possible important modifications, of the condi- 
tions (|61ti64l) to future work and focus in the following 
discussion on well-posedness of the constraint subsector. 
In the following sections conditions (|6HI64D are consid- 
ered only in a spherical reduction of the system. The 
constraint absorption properties of conditions (1601) are 
studied in our numerical tests. 



III. ANALYTICAL SETUP 

In this section we discuss the strategy adopted to prove 
well-posedness, namely the frozen coefficient approxima- 
tion. This simplification is necessary since one can show 
that our system is not symmetric hyperbolic and our BCs 
are not maximally dissipative. We also outline the cas- 
cade method which can be used for general proofs, pro- 
vided that the equations of motion of the system have a 
special structure. 

A. Frozen coefficient approximation 

Once the BCs are specified, one should determine 
whether or not the resulting IBVP from (I28H30|) with (|BD} 
164]) is well-posed. This can be done by using the 
frozen coefficient technique, where one considers a high- 
frequency perturbation of a smooth background solution. 
This regime is the relevant limit for analyzing the con- 
tinuous dependence of the solution on the initial data. 
By considering such a perturbation one can detect the 
appearance of high frequency modes with exponential 
growth. Therefore, if the IBVP is well-posed in the frozen 
coefficient approximation, it is expected that the prob- 
lem is well-posed in the non-linear case. In this limit, the 
coefficients of the equations of motion and the bound- 
ary operators can be frozen to their value at an arbi- 
trary point. So, the problem is simplified to a linear, 
constant coefficient problem on the half-space which can 
be solved explicitly by using a Fourier-Laplace transfor- 
mation 0, [431 - This method yields a simple algebraic 
condition (see appendix EJ which is necessary for the 
well-posedness of the IBVP. 



Following [23] , we perform a coordinate transformation 
which leaves E t invariant, such that one can rewrite the 
background metric (|53p at the point p in the form 



ds 2 \ 



dr+[Pdt + dx) +Ay z + dz\ (65) 



where P corresponds to the normal component of the 
shift vector at p. According with this approximation, 
one can assume that the BC is a plane. Therefore, the 
spatial manifold becomes £ = {(x,y,z) g R 3 : x > 0}. 
We denote the flat spatial metric at p by rjij. Regard- 
ing the above metric, the time-like and outgoing normal 
vectors (n°, s a ) are 



h a d a = d t -pd x , s a d a = -pd x 



(66) 



Besides, by using and choosing 1+log slicing hl — 
2/ a and fixing fig = 1 in the shift condition one can 
rewrite the equations of motion (|28H30I) in the frozen co- 
efficient approximation at a point p in the form 



($-2#d l )a = 0, 
(dl - d l d t ) p l = Q - O l d a 



(9 2 - 3%) 7ij = -(I -2e x )r, M d i d jlu 
+ 2(1 — e a )didja . 



(67) 



(68) 



(69) 



Note that with the additional choice of asymptotically 
harmonic shift (e a , e x ) = (1, 1/2) the resulting IBVP for 
the above system with boundary conditions 



(|60H64|) has 

a cascade property; the gauge sector (a, Pi) is coupled 
with the metric only through the BCs. One can analyze 
the resulting IBVP for the gauge sector and then use it 
as a source in the remaining system. Nevertheless, with 
the standard choices (e a ,e x ) = (0,0), all the variables 
are coupled to each other in the bulk as well as at the 
boundary. In this case, one should analyze the full system 
simultaneously. We will present the analytical results for 
arbitrary value of these parameters in [48j . 



B. 2+1 decomposition 

To rewrite the above system with BCs (|60M64p as a set 
of cascade of wave problems, we perform a 2 + 1 decom- 
position against the spatial unit vector s 1 = —e x . Thus, 
the lapse satisfies 



df - 2 p d t d x ~(2-p 2 )d 2 x -2 d A d A 

d t - (V2 + °P) d x ~\ L+1 a = h a . 



a = 0, (70) 
(71) 
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The wave problem for f3 s and (3 A is obtained once we 
project the system 



&t - 2f3d t d x - (1 - - d c d c A 



9( 



(l + /3j d x \ Pi = h 



-\d d ia , (72) 
(73) 



along s l or along the transverse directions, respectively. 
Note that the lapse is only coupled with the equation ([72"|) 
in the bulk. One could naively think in to analyze the 
lapse subsystem independently and use it as a source 
in ([72]). Therefore, the resulting global estimate for the 
gauge sector will contain more derivatives of the lapse 
than of the shift. In general, it can spoil the estimate 
when one considers lower order terms which appear in 
the nonlinear case. In order to prevent it one should 
consider the wave problems for the lapse and the shift 
simultaneously. 

The system for the 7^ is 



d*-2j3d t d x 



di 



(l + P) d x 



(1 - f)dl 

L+l 

„,TF j 
lAB = 



d c d, 



•c 



t2I = o, 



, TF 

n AB 



(74) 
(75) 



According to the results presented in the appendix |Bl 
it is straightforward to prove this problem is well-posed. 
Since the above subsystem is decoupled completely from 
the rest of the system, it can be considered as a given 
function in the remaining wave problems. On the other 
hand, by introducing the trace variable A and considering 
the equation pip as a definition of the constraint O, we 
obtain 



df - 2/3 dtd x - (1-/3 



1 + )d m 



2 )dl 

L 



d A d A 



A = 0. 



(d t -/3d x )(a-A) 



2d x p s + 2d A (3 A 



= 0. 



(76) 



(77) 



Since this system is also decoupled from the rest of the 
metric sector, one can analyze the resulting IB VP and af- 
ter that, use it as a given source for the other problems. 
Finally, the remaining equation of motions are again ob- 
tained through the projection of the wave equation 



d?-2pdtd x 



d c dr 



Us = , (78) 



along s l or along the transverse directions respectively. 
By virtue of Eqn. (02j the BC for j ss is 



di 



1 + /3 a 



{d t -pd x )p s -d A lsA 



d x (a - 7ss + A/6) 



(79) 



and finally, the BC for j sA is 

L r 



d, 



(l + /3) 



{d t -pd x )p A -d B 1 H + 



d xlsA + d A (a- A/3 + lss /2) = . (80) 



This subsector does not have the cascade property. The 
equation of motion for j ss and j sA are decoupled, but 
their BCs are coupled to each other. Therefore, one 
should consider the mutually coupled wave problems ([751 
[50)1 simultaneously. 

In the following section, we consider the case with 
trivial initial data. Notice that this is not a real re- 
striction. One can always treat the case of general ini- 
tial data by considering, for instance, the transformation 
u(t, x) = u(t, x) — g(t) f(x), where g(t) is a smooth func- 
tion with compact support such that g(0) — 1. Therefore, 
u(t,x) satisfies the same IBVP as u(t,x) with modified 
sources and trivial initial data [2|. 



IV. WELL-POSEDNESS RESULTS 

This section contains our analysis of well-posedness for 
different BCs. As we have mentioned before, we consider 
the IBVP for the constraint subsystem on the manifold 
M = [0, T] x S and we just analyze the correspond- 
ing IBVP for the dynamical Z4c variables on a spheri- 
cal scenario. To do this, we explicitly solve the bound- 
ary problem using the Laplace-Fourier transformation. 
Kreiss presented in 6] sufficient conditions for the well- 
posedness in the frozen coefficient approximation. The 
key result in Q is the construction of a smooth sym- 
metrizer for the problem for which well-posedness can be 
shown using an energy estimate in the frequency domain. 
We summarize the principal ideas of the Kreiss theory in 
the appendix [A"l 



A. Constraint subsystem 

Consider the IBVP for O with first order BCs (L=0) 
and, for a moment, let us assume inhomogeneous BCs, 
i.e. = q, where q is a given boundary data. One can 
show that the equation of motion (|33[) and the boundary 
for this variable can be rewritten in the form 



(s 2 +Lu 2 )~2psd x -{1- f3 2 )8 2 x 6 = 0, Q = q 



where O and q denote the Laplace-Fourier transforma- 
tion of O and q with respect to the directions t and x A 

respectively and u> - 



oj 2 + uj 2 . Following jia |23J we 



rewrite the above system in the form (IB6tiB7|) by intro- 
ducing the variable 



dq = - (d x e + 7 z u p s e 



(81) 



with W 



(e,£e) T , L(a,u) = (1,0) 



and 



M(s,u) = K 



-j 2 /3s' 

„4 \2 



7 4 A^ - 7 2 /3s' 



(82) 
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where we have defined 7 = 1/y 1 — /3 2 , k = -\/|s| 2 + w 2 , 
s' = s/k, l>' = uj/k and A 2 = s' 2 + 7 ~ 2 w' 2 . The corre- 
sponding eigenvalues and eigenvectors of M = M(s,w) 
are 

r± = -k 7 V/3tA), e± = (1,± 7 2 A) t . (83) 

Using this, it can be shown that the L2 solution of the 
svstem (|B7p is given by 



W(s, x, lu) — a e exp(r x) 



(84) 



where a is a complex integration constant which is de- 
termined by introducing (|84[) into the boundary, i. e. this 
constant satisfies a — q. Therefore, it follows immedi- 
ately that 



\W(a,0,u)\<C\q\ 



(85) 



where C is a strictly positive constant. Provided that the 
eigenvectors in the solution (|8~4"]) are normalized in such 
a way that they remain finite as U) — >■ 0, as uj — > ±00 and 
as |s| — > 00 then we conclude that the resulting IB VP 
for (|3"3"|) with the first order boundary condition ([50)1 
is well-posed for trivial initial data. By inverting the 
Laplace transformation and using the Parseval's identity, 
we obtain an hi estimate of the form 



\\W(;t)\\ldt<C T I hWl^dt., 



(86) 



in the interval < t < T for some strictly positive con- 
stant Ct > 0. Note that if the constraint O is satisfied 
initially and we consider trivial boundary data q = 0, 
the above inequality implies that the constraint is satis- 
fied everywhere and at each time. An equivalent estimate 
for Zi holds. 

I 



We generalize our previous analysis to consider BCs 
which depend on second or higher order derivatives of 
the constraints. Recently, it has been shown that those 
conditions reduce the amount of spurious reflections at 
the boundary [HIH- In fact, the BCs (|BDJ| are Dirich- 
let conditions for the constraint subsystem (f3"3"|) . which 
means that the constraint violations that leave the bulk 
are reflected at the boundary. Therefore, Let us consider 
the wave problem for the O with high order BCs [L > 1). 
Note that according to the appendix iBl it is possible to 
rewrite the boundary matrix L = L(s,ui) in the form 



A7 2 



(87) 



where a± = s'±A. The L2 solution of this system is given 
by ((84)) but now the integration constant a satisfies 



ah cr 



(88) 



According to appendix [B] one can show that there is 
a strictly positive constant S > such that |a + | L = 
\s' + A| > S. Therefore, it follows that the system sat- 



isfies the estimate ([85]) . We conclude that the solution 
of the system (|33[) with high order constraint preserving 
BCs dSOJ is well-posed and it satisfies ([86]). The IB VP 
for Zi with these BCs can be treated in a similar manner. 
B. Gauge subsystem 

Consider now the spherical reduction of the wave prob- 
lems (|70H73[) for the lapse and the shift with first order 
BCs, i.e. conditions ([7T|) and (1731 with L = 0. By intro- 
ducing the first order reduction variables 



Da=- [d x + 1 z a (3s)a 



(89) 
(90) 



it can be possible to rewrite those wave problems in the 
form (IB7I) with 



K V 



W 



L(s) = 



( a \ 
Da 

( V2s 

\ 



, M(s) = K 
~T a 2 



1 

-llfrs' 

-7 2 /3s' 1 

V-2s' 2 /3 7 2 7 4 s '(2 + /3 2 ) 7 2 7 2 /2 7V 2 - 7 2 /3s'/ 

9s J K\(l-p)h s 



-Y 4 S' 2 
la °a 



9 = 



r 



Here we have defined s' = s/k, k — \s\ and 



The I/2 solution of the above system is given by 



2-(3 2 



1-/3 2 



(91) 



W(s,x,u) 



i=l 



a i e l exp(?v x 



(92) 
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where oi are the complex integration constants and r~ 
and e~ the negative eigenvalues of M and e~~ the corre- 
sponding eigenvectors. By replacing the solution into the 
boundary condition, it is possible to show that <7 2 : satisfy 



(93) 



9a 



(1-/3) (2 + V2 + ^(l + V2) 



2 s' 



8 s' 



(94) 



By the same argument as above, it is straightforward by 
using the triangle inequality that the system is boundary 
stable and therefore an equivalent estimate as (l86t holds 
for the gauge subsystem with first order BCs. 

Next, let us consider high order BCs for this subsys- 
tem. To give an estimate for the solution of the gauge 



subsystem with high order conditions we rewrite them 
in an algebraic form. Therefore, by using the procedure 
presented in appendix [B] and by virtue of the equations 
of motion (I67H68P , one can rewrite the boundaries (|72j|73[) 
with L = in the form CW = AW where 



A = 



( V2s' 




-7« 
V2s' 




7 



\ 



2 



\2s' 2 M 4 a s'(2 + /3 2 ) 7 2 /2 -7 2 * 2 sf ) 



We have defined the boundary operator by C — 
£i) T with 



£ M = (fx - /§) s' 



d x 



By iteration, explicitly one obtains C L+1 W 
Here we have defined 



(95) 
A L+1 W. 



a l+i = 



a L+1 
-V2s^ a a L+1 



-a^/i^s'jl) 



-F0) li «+ ++1 /7 2 G0) 7 2 a L +*/(s' 7 2 ) 
V H{°P)s' 1 ia L + +1 -J(« 7 ^t +1 



-s 7 






b L + +1 
H L + +1 






-^ +1 /( S '7 2 ) 



where a + = 2 \[2 s' and b + = 2 s' are the eigenvalues 
of the matrix A = A(s') and F((3), • • • , J($) are certain 



shorthands for combinations of /3 only. Thus, the high 
order BCs for the gauge subsystem can be rewritten in 
the form (|B7[) with the boundary matrix given by 



L(s,u)) = 



1 



a L + +1 -a L + 
^(/3)7^^ +1 /7 2 G(/3) 7 > i +VGs'7 2 ) 



-V(V2 S ' 7 2 ) 



lL+1 



-6 





'V( S '7 2 ) 



(96) 



and the given data 



((V2-P) L+1 h c 
{ (l-/3) i+1 ^ 



(97) 



The L2 solution of the system is given by f|92|) . Never- 
theless, the complex integration constants now satisfy 



fjo 



2a L + 1 ' 

>(/§) + V2 G{P)lt 2(l + V2)(l-/3) 



4 6 i+1 7 2 



8 (2 - V2/3) a 



L + l 



9s 



2b L + 1 ■ 



(98) 



(99) 



Since Re(s') > 0, 



-,L+l 



and b L+1 are proportional to s' 



and the remain coefficient are constants, it follows that 
the system is boundary stable and therefore well-posed. 



Similar arguments apply to the wave problem of the met- 
ric components. 



V. NUMERICAL APPLICATIONS 

The necessity of CPBCs for the Z4c system is not only 
motivated by the fundamental requirement of having a 
mathematically well-posed system, but also by the nu- 
merical evidence of artifacts related to the implementa- 
tion of inadequate BCs. The property of full propagation 
of the constraints is both a strength and a weakness of 
the Z4 evolution system. On one hand, by a compari- 
son with BSSN it was shown in [37] that this property 
reduces constraint violations on the grid. On the other 
it makes the BCs a more important issue, because if the 
numerical boundary condition introduces large constraint 
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FIG. 1. Radial oscillations of the central rest-mass density 
of a compact star in time. When Sommerfeld conditions are 
used a perturbation from the boundary (at r out = 20) hits 
the center of the star (at t = 20) and further perturbs it. 



violations, perhaps as spurious reflections, then the vio- 
lation may propagate inside the domain and swamp the 
numerical solution. 

As an example of numerical artifact, pointed out but 
only briefly discussed in [iST] . we show in Fig. [TJ the time 
evolution of the central rest-mass density of a equilib- 
rium model of spherical compact star obtained with the 
Z4c formulation (coupled to general relativistic hydrody- 
namics equations) and two different BCs: Sommerfeld 
and CPBCs [5(|. The central density is expected to re- 
main constant in time but the truncation error of the 
numerical scheme causes small oscillations around the 
initial value that converge away with resolution. This is 
clearly visible when CPBCs are used. The frequency of 
these oscillations corresponds to the proper radial mode 
of the star. When Sommerfeld conditions are used how- 
ever the constraint violation at the boundary propagates 
into the domain and perturbs the star as soon as it be- 
comes causally connected with the outer boundary. The 
perturbation from the boundary alters the mean value of 
the central density. The "boundary" perturbation does 
not converge (or converges much slower than the inte- 
rior, see Sec. IV Bp so becomes the dominant error when 
the grid is refined. At later times the oscillations are 
damped by the hydrodynamical interaction between the 
fluid and artificial vacuum, but the mean value of the 
central density is forever modified from the initial one. 
Such boundary artifacts can thus dramatically move the 
solution far away from the initial configuration in phase 
space, despite the fact that at late times the constraint 
violation is very small. Preliminary 3D simulations with 
the Jena BAM code [13] of single neutron stars, single 
puncture and binary black holes showed features very 
similar to the spherical results and in some cases insta- 
bilities triggered by the boundary. 



In the following sections we discuss numerical results in 
spherical symmetry focusing only on the boundary con- 
ditions ((60)) . Since the practical implementation of BCs 
in a code is also an issue, in Sec. IV Al we summarize the 
method we used as well as other standard approaches. 
We use the code of [371 ]. For more information on our 
numerical method, spherical reduction of the equations 
please refer to appendix A of that reference. We per- 
form several tests, in each case with Sommerfeld and con- 
straint preserving conditions. To examine stability and 
the effect of the BCs we consider simulations with a very 
close outer boundary (r out ~ 20 M) and compare the re- 
sults with a reference simulation [231 l5l| , in which the 
outer boundary is placed far away {r' out ~ 1000 M) from 
the origin. Since the boundary of the reference solution 
is causally disconnected on the time scale of the simula- 
tions with closer boundary, the BCs in the reference run 
have no importance. Results are presented for moderate 
resolution, about Ar ~ 0.12; but higher resolution runs 
as well as very long-term (hundreds of thousands of cross- 
ing times) simulations were performed showing the same 
behavior and no numerical instabilities. To monitor the 
global constraint violation we define the quantity: 



C = \JH 2 + M i M i + 6 2 + Z l Z, , 



(100) 



and we will refer to it as the constraint monitor. We will 
often make use of 2-norms of quantities: 



l|C(-,t)|| a 



drr 2 C{r,tf 



(101) 



where in practical computations the integral is performed 
on the grid by the trapezium rule. For a fair comparison 
with the "near-boundary" solution, the norm of the ref- 
erence solution is taken only on the domain, r £ (0, r out ). 

Since most of the analytical results were obtained with 
the new asymptotically harmonic shift condition, we ex- 
amine this gauge as well as the standard puncture gauge. 
In all cases we have found comparable results (see e.g. 
Fig. [4]) and therefore we will focus primarily on the stan- 
dard puncture gauge. We therefore aim to give at least 
some numerical evidence of well-posedness in those cases 
where we are unable to demonstrate strong mathematical 
results. 

A brief description of the tests performed and their 
aim follows below. 

Perturbed flat spacetime. Evolution of constraint vi- 
olating initial data on flat space. Here we focus on con- 
vergence and constraint absorption. We find near-perfect 
constraint transmission of the constraints when using the 
second order CPBCs. We devote the most attention to 
this test because the effect of the BCs are clearest in the 
absence of other sources of error. 

Star spacetime. Evolution of a stable compact star. 
In the Sommerfeld case, non convergent reflections from 
the boundary effect the dynamics of the star. The ab- 
sorbing CPBCs completely solve this problem. 
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Black hole spacetime. Evolution of black hole initial 
data. The robustness and performance of CPBCs are 
tested against black hole spacetimes with different ini- 
tial data and gauges. In particular we evolve a single 
puncture and Schwarzschild with a Kerr-Schild slicing. 

We use geometrical units c = G = 1 everywhere, in 
case of matter spacetime dimensionless units are adopted 
setting Mq = 1, in case of black hole spacetime Mbh = 
1, while in case of fiat spacetime a mass scale remains 
arbitrary. 



A. Numerical implementation of boundary 
conditions 

The literature contains many suggestions for the im- 
plementation of BCs, of which we highlight a small subset 
here. 

Populate ghostzones. One example is the recipe 
of [52| . under which numerical stability of the shifted 
wave equation was proven. The idea is to write the BCs 
on the first order in time second order in space char- 
acteristic variables (which only implicitly contain time 
derivatives) and populate ghostzones so that the desired 
continuum boundary condition is satisfied. Ghostzones 
not determined by the BCs are simply populated by ex- 
trapolation. Since this method relies in an essential way 
on altering spatial derivatives locally, it is not clear how 
to apply the approach if one is computing spatial deriva- 
tives pseudo-spectrally. Furthermore the recipe may not 
give a unique prescription when one is given a system of 
equations. 

Summation by parts. As we have mentioned before, 
with the summation by parts schemes and penalty tech- 
niques a quantity that mirrors the continuum energy for 
a symmetric hyperbolic system is constructed on the dis- 
crete system. Since we do not rely on an energy (sym- 
metric hyperbolicity) method in our well-posedness anal- 
ysis we are not able to construct a summation by parts 
finite difference scheme that guarantees stability. Analy- 
sis and applications in numerical relativity can be found 
in [ilH^I. 

BCs as time derivatives of evolved fields. The remain- 
ing approach we consider assumes that one may take 
spatial derivatives everywhere in a spatial slice. Inside 
a pseudo-spectral method spatial derivatives are natu- 
rally defined everywhere. On the other hand when ap- 
proximating derivatives by finite differences, one may ei- 
ther use extrapolation to populate ghostzones, or take 
lop/one-sided differences near the boundary. The two 
methods are equivalent. To express the BCs on time 
derivatives of the metric one starts by substituting the 
definitions of the time-reduction variables (for example 
<d,Zi and Kij) into the BCs. Then, if higher than first 
order time derivatives are still required, one may define 
a set of auxiliary reduction variables. One may use the 
equations of motion of the auxiliary variables in combi- 
nation with the boundary conditions to eliminate spatial 






50 100 150 200 250 300 350 400 
Time 

FIG. 2. Constraint violation in flat spacetime test. The 
2-norm of the constraint monitor is showed in time for differ- 
ent BCs implemented. The same quantity for the reference 
simulation is showed. 



derivatives of auxiliary variables so that they may be con- 
fined to the boundary [45], [58| . This approach was tested 
in the Caltech-Cornell SpEC code (24|, but it is currently 
not used for binary black hole simulations. Here we per- 
form evolutions only up to second order, and so do not 
need to use auxiliary variables in the boundary. 

Our numerical tests are performed with a spherical re- 
duction of the Z4c system. The numerical implementa- 
tion of the boundary conditions in spherical symmetry is 
described in appendix [Cj 



B. Perturbed flat spacetime 

In this test an initial constraint violating Gaussian per- 
turbation is prescribed in variable x- During the evo- 
lution it propagates, reaches the boundary and if com- 
pletely absorbed, the system relaxes to the Minkowski 
solution. This does not happen in practice because re- 
flections are always present. Here we investigate the mag- 
nitude of these reflections and compare CPBCs of 1st and 
2nd order with standard Sommerfeld conditions. 

Figure [5] shows the 2-norm of the constraint monitor 
in time, results from the reference simulation are also re- 
ported. All the data agree up to around t — 20, i.e. when 
the solution is traveling through the grid. After that time 
the following happens: (i) the constraint violation re- 
mains almost constant for 1st order CPBCs (green dotted 
line) and eventually causes the simulation to crash; (ii) 
the constraint violation decreases for Sommerfeld (blue 
dashed line) initially not monotonically (notice the four 
plateaus) then, after t = 180, reaching a monotonic be- 
havior; (iii) the constraint violation for 2nd order CP- 
BCs (continuous red line) is smaller than Sommerfeld, 
plateaus are less clear but a monotonic behavior is also 
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reached after t = 180; (iv) the reference solutions (black 
dotted line) agree among themselves (remember that on 
the plotted timescale the boundary of the reference solu- 
tion is disconnected from the spatial domain under con- 
sideration) and decrease monotonically but in absolute 
value less than 2nd order CPBCs. 

The interpretation of these observation is quite obvi- 
ous, at least for points (i)-(iii): 1st order CPBCs sim- 
ply do not absorb the perturbation, which is entirely re- 
flected back; Sommerfeld BCs are affected by partial re- 
flections and at each crossing time a smaller portion of 
the wave comes back into the domain until an almost 
complete absorption; 2nd order CPBCs immediately ab- 
sorb the largest portion of the outgoing wave. These 
statements are quantified below. 

Let us briefly discuss point (iv). A close look at the 
evolution of the constraint monitor in space shows a sys- 
tematic drift from zero due to a back-scattering effect 
of the perturbation leaving the grid. This is responsible 
for the larger values of the constraint violation. Such 
an effect can not be simulated in the closer domain case 
because the CPBCs cut completely the incoming modes 
of the solution, and, effectively, two different numerical 
spacetimes are simulated in the two cases. 

An important point is to quantify the absorption prop- 
erties of the BCs under consideration. To this end we 
consider the characteristic fields associated to the given 
by: 



d e±e.. 



(102) 



They can be regarded as incoming and outcoming modes 
of the solution. Thus we define an experimental reflec- 
tion coefficient (inspired by [46l l49l|) defined as the ratio 
of the Fourier modes of the characteristic fields at the 
boundary: 



R 



_ \Ue(k)\ 

" \u£(k)\' 



(103) 



Figure [3] shows that R ~ 1 for 1st order CPBCs, while 
the behavior of 2nd order CPBCs is qualitatively similar 
to Sommerfeld. Since they do not absorb the constraint 
violation, in what follows, we discard the 1st order CP- 
BCs. 

The results presented so far refer to the puncture 
gauge. As stated at the beginning of the section, basically 
no significant differences are found when the asymptoti- 
cally harmonic shift is employed. Figure S] shows clearly 
this fact, we do not further comment on the asymptoti- 
cally harmonic gauge until Sec. IV Dl 

Finally, we present convergence results. In figure[5]the 
experimental self-convergence factor is plotted in time. 
While the physical solution is traveling on grid, t < 20, 
the scheme is fourth order convergent as expected. After- 
wards the numerical solution consists of the boundary re- 
flections only and, while in case of Sommerfeld reflections 
are first order accurate, in case of CPBCs they maintain 
fourth order convergence up to t ~ 100. For later times 
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FIG. 3. Experimental reflection coefficient in flat space- 
time test. The experimental reflection coefficient, defined in 
Eq. (|103|) . is plotted versus the wave number for different BCs 
implemented. 
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FIG. 4. Constraint violation in flat spacetime test. Compar- 
ison of 2nd order CPBCs with puncture and asymptotically 
harmonic shift. 



(not shown in the plot) the absolute value of the solution 
(and of the reflections) is so small that only noise is seen. 
We remark that the order of extrapolation used to fill 
the ghosts points and the finite difference operators in 
the two approaches are the same, so the differences are 
really due only to BCs. In order to obtain these conver- 
gence results, a non-staggered grid must be used because 
the staggered grid converges to the continuum domain at 
only first order in the grid spacing. In the following ap- 
plications we use staggered grids, since this is commonly 
done in 3D codes. 
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FIG. 5. Convergence factor in flat spacetime test. The self- 
convergence factor is compute from the 2-norms of three sim- 
ulations at different resolutions and showed in time for 2nd 
order CPBCs and Sommerfeld BCs. 



FIG. 6. Constraint violation in star spacetime test. Upper 
panel: The 2-norm of the constraint monitor is plotted in time 
for different BCs implemented. The same quantity for the ref- 
erence simulation is shown (black dotted line). Bottom panel: 
The 2-distance of the constraint monitor with the reference 
simulation is showed in time for different BCs implemented. 



C. Star spacetime 



D. Black hole spacetime 



In this test we evolve stable spherical star initial data 
(see [13] for detail). At the beginning of the section we 
already discussed one of the main drawback of the use 
of Sommerfeld BCs with Z4c. As it is evident from Fig- 
ure m the use of CPBCs is not optional but absolutely 
necessary to obtain reliable results. In a more compli- 
cated/dynamical scenario in fact such artifacts could be 
hidden or erroneously interpreted as real physics. 

In this paper we are presenting results obtained with- 
out the Z4 constraint damping scheme [36| . One may 
however consider using the constraint damping scheme 
with a sufficiently large computational domain to sup- 
press perturbations from the boundary. In our experi- 
ence this is not only a inefficient cure (especially in 3D 
simulations) but an ineffective one. Our simulations in- 
dicate that the required damping coefficients are quite 
large, possibly because the perturbation from the bound- 
ary is typically not a high-frequency perturbation, and 
the damping scheme is most effective on high frequency 
perturbations. Constraint damping is analytically under- 
stood in the linear regime and high frequency approxima- 
tion, in a more general situation the indiscriminate use of 
constraint damping may lead to undesirable effects (e.g. 
qualitatively similar to large artificial dissipation). 

In the upper panel of Figure [5] we show the 2-norm of 
the constraint monitor for Sommerfeld, 2nd order CPBCs 
and the reference solution. It is evident that CPBCs are 
closer to the reference solution. In the bottom panel the 
distance from the reference solution is plotted showing 
that CPBCs lead to an improvement of approximately 2 
to 4 orders of magnitude. 



In this test wc consider different black hole space- 
times. A spherical black hole was evolved with punc- 
ture and Kerr-Schild initial data; evolutions were per- 
formed with both the gamma driver and asymptotically 
harmonic shift. 

We focus first on the evolution of a single puncture. 
Figure [7] shows that the behavior of the constraint mon- 
itor, computed outside the apparent horizon, is analo- 
gous to what we found in the flat and star spacetime 
tests. As the solution asymptotes to the stationary trum- 
pet slice [59l | we find more constraint violation than in 
the evolution of the star. The resolution employed in 
the simulations for the figure is quite moderate and the 
outer boundary very close, even compared with a 3D 
code. Nonetheless the CPBCs perform quite well and 
significantly improve the numerical solution with respect 
Sommerfeld. No big differences were found when adopt- 
ing either the gamma driver shift or the asymptotically 
harmonic one. As shown in Figure [5] the puncture shift 
(red solid line) does even better than the asymptotically 
harmonic (red dashed line) after t ~ 60. To the best 
knowledge of the authors this is the first time that the 
asymptotically harmonic shift has been used in the evo- 
lution of puncture data. 

To further assess the robustness of our BCs we evolve 
the spherical black hole with Kerr-Schild initial data: 



ds* = - I 1 - — I dt 2 + —dtdr 



+ r 2 dn 2 . 



1 + — J dr 2 
(104) 



The excision surface is at r = 1.9M and simple extrapo- 
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FIG. 7. Constraint violation in black hole spacetime test. 
Upper panel: The 2-norm of the constraint monitor is showed 
in time for different BCs implemented. The same quantity for 
the reference simulation is showed (black dotted line). Bot- 
tom panel: The 2-distance of the constraint monitor with the 
reference simulation is showed in time for different BCs im- 
plemented. 



lation is used at the inner boundary. We stress that this 
test is initially more demanding for the gamma driver 
shift. Since the line-element is not written in a manifestly 
conformally flat form, the shift immediately evolves ev- 
erywhere in space, which is not the case with the gamma 
driver shift and puncture data. Figure [8] shows that the 
results for the constraint monitor obtained with our CP- 
BCs are comparable to the puncture case with both the 
puncture shift (orange thick solid line) and the asymp- 
totically harmonic (orange thick dashed line). In this 
case the latter perform better. Moreover, at this resolu- 
tion, the bigger spurious reflections produced by Som- 
merfeld conditions combined with the close boundary 
used (r out = 20) cause the simulations to crash. To 
achieve stable evolutions with the Sommerfeld conditions 
we find it necessary to use a higher resolution and a more 
distant outer boundary. 



VI. CONCLUSION 

For numerical applications of free-evolution schemes 
in general relativity there is a strong motivation to con- 
struct constraint preserving boundary conditions. With- 
out such conditions, constraint violations may appear at 
the boundary of the numerical domain and swamp the 
numerical solution. In the worst case such errors could 
be interpreted as real physics. In this study we have 
constructed constraint preserving conditions for the Z4c 
formulation of the Einstein equations with variations of 
the popular puncture gauge. 

We demonstrated well-posedness of the resulting ini- 



FIG. 8. Constraint violation in black hole spacetime test, 
comparison of different initial data and gauges. The 2-norm 
of the constraint monitor is showed in time for 2nd order CP- 
BCs. The evolutions refer to puncture initial data showed 
with puncture shift (red solid line) and asymptotically har- 
monic shift (red dashed line) and to Kerr-Schild initial data 
(evolved with excision) with puncture shift (thick orange solid 
line) and asymptotically harmonic shift (thick orange dashed 
line). 



tial boundary value problem for the constraint subsys- 
tem on a four dimensional compact manifold in the high- 
frequency approximation. Since we are only interested 
in the constraint absorption properties of our boundary 
conditions we just analyzed the initial boundary value 
problem of the a spherical reduction of the Z4c system 
for a special choice of free parameters of the gauge con- 
dition. One may be able to expand our calculations for 
the asymptotically harmonic shift condition by verifying 
the existence of a suitable symmetrizer in a neighbor- 
hood of the point p about which we perturb to reach the 
high-frequency limit. Unfortunately there is not a clear 
method for extending our calculations beyond the high- 
frequency limit with the standard gauge choice. Even 
in the high-frequency approximation we find that the 
Laplace-Fourier method becomes cumbersome. The key 
problem is that the cascade approach of Kreiss and Wini- 
cour fails with the puncture gauge. 

In order to build a body of evidence for the well- 
posedness of the initial-boundary value problem with 
the standard puncture gauge we therefore performed nu- 
merical evolutions of various spherical initial data sets. 
Roughly speaking our approach for the implementation 
of the boundary conditions is to rewrite them as closely as 
possible to Sommerfeld conditions most commonly used 
in BSSN evolutions. We will report our method in detail 
elsewhere. Since the underlying formulation is not sym- 
metric hyperbolic we are not able to make a summation- 
by-parts approach to the implementation. Therefore nu- 
merical stability is most straightforwardly established by 
studying toy problems mathematically and performing 
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existence numerical tests on the numerical system. We 
compared the behavior of the system with the standard 
puncture gauge and the asymptotically harmonic gauge. 
We find very similar features in all tests. We demon- 
strated that with our approach to the numerical imple- 
mentation we can achieve clean pointwise convergence 
even in reflected constraint violation. We also found, in 
agreement with previous studies, that high-order bound- 
ary conditions are able to absorb outgoing constraint vi- 
olation much more effectively than first order conditions. 

There are two obvious places in which we would like to 
strengthen our results. Firstly it is desirable to extend 
our well-posedness results, especially in the case of the 
standard puncture gauge. However in order to do this a 
different mathematical approach will probably be neces- 
sary. Secondly, that the numerical tests were performed 
in spherical symmetry is an obvious drawback which we 
aim to address shortly. Preliminary tests in 3D with the 
BAM code indicate that additional tangential terms in 
the boundary conditions, which we have discarded in this 
paper, are required in order to get a stable evolution. 
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Appendix A: Well-Posed problems 

The well-posedness of the IBVP is the requirement 
that for given initial and boundary data an unique so- 
lution should exist and it should depend continuously on 
the data (see e.g. [6(| )■ In this section we present a short 
review of the theory for first order systems to prove well- 
posedness This theory developed by Kreiss [y] gives us 
necessary and sufficient conditions for well-posedness of 
the IBVP for strictly hyperbolic systems. The theory was 
extended to hyperbolic systems of constant multiplicity 
by Agranovich [7J • The discussion presented here follows 
closely that of in @, [l§ HI- 

Consider a strongly hyperbolic system of equations of 
motion 

n 

dtu = B d x m + C A d A u , (Al) 

A=2 

with constant coefficients on the half-space t > 0, x 1 > 
and — oo < x 2 , . . . , x n < oo. Here u is a n-dimensional 



vector and B and C l are n x n constant matrices. By 
assuming that B is non-singular, it can be rewritten as 

B = ( ~f A °„ ) , (A2) 

with A 1 and A 11 real and positive definite diagonal ma- 
trices of order to and n — m respectively. We imposed m 
BCs at x 1 = in the form 

Lu'fax)}^ =g{t,x), (A3) 

where L is a to x to constant matrix and g — 
g(t 7 x 2 , . . . , x n ) is a given boundary data vector. For sim- 
plicity, we consider trivial initial data u(0, x) = 0. In the 
following we solve the IBVP (|A1HA3|) by performing a 
Laplace-Fourier transformation with respect to the di- 
rections t and x A tangential to the boundary x 1 = 0. 

Let u — u(s. x 1 , oj a ) denote the Laplace- Fourier trans- 
formation of u(t,x). Then, u satisfies the ordinary dif- 
ferential system 

d x u = M(s,u>) u , oni 6 (0,oo), (A4) 
Lu 1 = g, at x = , (A5) 

where g denotes the Laplace- Fourier transformation of g 
and 

M(a,u) = B- 1 (sI nxn + iu; A C A ). (A6) 

If Xi and e.j(s, u>) are the corresponding eigenvalues and 
eigenvectors of M then the general solution (functions 
which are quadratically integrable) of (|A4[) is given by 

m 

u = ^ff.e,(s,w) exp(A i x 1 ) , (A7) 

i=l 

where <7j are complex integration constants [6l| . These 
constants are determined by the boundary conditions. 
By substituting l|A7[) into the expression (|A5[) we obtain 
a system of m linear equations for the unknown cr,. This 
system can be written in the form 

B(s,u)a = g, (A8) 

where D(s,w) is a to x m matrix. Let us consider for a 
moment homogeneous BCs g = and suppose that there 
is s with Re(s) = rj > such that DetO = 0. It means 
that (|A8[) has a non-trivial solution and therefore, the 
solution of the IBVP (|AHA3|> is 

u(t, x % ) =u(x 1 )exp(st + icu A x A ) . (A9) 

This implies, by homogeneity of the system (|A4IIA5|) . that 

=u(9x 1 )exp((st + i(u) A x A ) . (AlO) 

is also a solution for any constant C > 0. By increasing 
that constant arbitrarily one can find a solution which 
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grows exponentially. Therefore, the IBVP is not well- 
posed. We conclude that the so-called determinant con- 
dition 



Dct © ^ . 



for 



ij>0, 



is a necessary condition for well-posedness. 

Next, let us consider the inhomogeneous BCs. Since 
the determinant condition is satisfied we can solve (|MJ 
for the integration constants. What remains to be shown 
is that solution (IA7|) can be bounded in terms of the data 
given at the boundary, 



|«(«,0,w)| <C\g(a,u)\ 



(All) 



where C > is an independent constant of s and uj. 
Using (jAllj) and by inverting the Laplace-Fourier trans- 
formation it is possible to show that the estimate 



T \\u(t 7 -)\\ 2 dt<6 f 
o Jo 



\g(t,-)\\ 2 \ xl=0 dt, (A12) 



holds, i.e. we can estimate the Li norm of the solution 
in terms of the Li norm of the given boundary data. 
Here the constant 6 is independent of the boundary data. 
Systems whose solution (IA 12|) satisfies this estimate are 
called boundary stable [2L [ a. |47l| . 

Kreiss has shown in [ft UM tnat if the system ([AT]) 
is strictly hyperbolic, the condition (IA11[) implies that 
there is a symmetrizer H — H(s / ,uj / ) such that 

I. H(s',uj') is a smooth bounded function that de- 
pends of (s',w'), 

II. there exists a constant S > such that 

HM + M* H> 5r)I, 

for all n > and all ui £ K, 

III. there are constants 62 > and C > such that 

(u,Hu) > 5 2 \u\ 2 - C\~g\ 2 , 

for all u satisfying the boundary condition (|A8|) . 



where s' — s/k, uj' = uj/k and k = \/|s| 2 + uj 2 and 

uj 2 = uj 2 +uj 2 . We have denoted the standard scalar 

y • 

product by (, • ,) and | • | its corresponding norm. This 
symmetrizer allows us obtain an estimate of the solution 
of the system for which we add a source term F(t,x) in 
the right of (|A1[) . In particular, according to [19| . we 
obtain an estimate of the form 



\\u(t,-)\\ 2 dt+ [ \\u(t,-)\\ 
Jo 



2 L=o dt < 



\F(t,-)\\ 2 dt + 



2 L=o dt 



(A13) 



where the constant S > is independent of the boundary 
data g or the source term F. 



Appendix B: Toy Model 



Consider the wave equation 

[dt-^&dt] U(t,x l ) 



F(t,x l 



(Bl) 



on the half-space t > 0, x > and y and z € (—00, 00) 
with trivial initial data. We impose BCs at a; = of the 
form 



[d a -fid x ] U(t,x l ) = h, 



(B2) 



where h are given boundary data and do is the time 
derivative along the coordinate time in the linear regime. 

According to [23], let us denote to U as the Laplace- 
Fourier transformation of U(t, x l ) with respect to the di- 
rections t, y and z tangential to the boundary then in 
the background (|6"5]) . U satisfies 



■ 2 /3 2 )d 2 x + 2p S d x -(s 2 + n 2 u J 2 ) 



(m 2 - /3 2 ) 



U = F; 



U 



(B3) 
(B4) 



where F and h denote the Laplace-Fourier transforma- 
tions of F and h respectively. In order to apply the the- 
ory presented in the appendix \^ one can rewrite the 
above system as a first order one by introducing the vari- 
able da m 



DU 



(B5) 



where 7^ = 1/y /1 2 — fi 2 . Therefore, the system (|B3|) can 
be rewritten in the form 



d x W = M(s,w)W + f , 
where we have defined 



(B6) 
(B7) 



W 



u 

DU 



K 



and 



M(s,oj) = k 



L(s,u) = Os',-7 2 ) , 



(B8) 
(B9) 



ifj, 

with A 2 = s' 2 + t,7 2 uj' 2 . The L2 solution of the homoge- 



■ - ■ +Jm ^ ■ ■ " - 

neous system (|B6|) is given by 



W(s, x, uj) 



(BIO) 



where r is the eigenvalue of M with Re(r ) < and 



e its corresponding eigenvector. Introducing (|B10|) into 
the boundary we have 



+ 7m 2 ^' 2 



(Bll) 
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According to [19j, 123J , one can show that there is a con- 
stant S > such that 



7m 2 



>S. 



(B12) 



Provided that the eigenvector e~ in the solution (IBlOp is 
normalized in a way that it remains finite as u> or \s\ goes 
to zero or infinity, there is a constant C > such that 



\W{a,0,u,)\ < C\g\, 



(B13) 



for all s € C with 77 > 0, and w e M. Therefore, the 
system is boundary stable. So, we should consider the 
existence of a symmetrizer = H(s',lu') and use it to 
get an energy estimate for the full problem. According 
to [23| , it can be shown that the following estimate 



DC 

V J (k(7| 2 + |^?7| 2 ) dx + (\nU\ 2 + \d x U\ 2 ) 







2 




1 


F 


dx + 


h 



x=0 



(B14) 



holds for some constant C > 0. Therefore, by inverting 
the Laplace-Fourier transformation and using the Par- 



seval's relation, one obtain the estimate (IA13|) for the 
solution in terms of the L2 norm of the boundary data. 

One can generalize the boundary condition (IB2|) to 
higher order BCs. It has been shown that such conditions 
reduce the amount of reflections at the boundary |24Ll49j. 
Thus, we impose high order BCs of the form 



[do 



fid x ] m+1 U(t,x*) 



(B15) 



with m > 1. Following [23J it is possible rewrite the 
previous conditions as 



£ m+1 U 



L+l 



(B16) 



where the linear operator C is defined by C — (fi — 0) s' — 
&x/kj 2 . Using the equation of motion (|B6|) . we can 
rewrite the above condition in algebraic form. Note that 



where the matrix A is given by 



.4 



H s 



k 2 I F 



-1u 



(B17) 



(B18) 



It has been shown in [23| that by iteration the boundary 
condition (IB16|) can be rewritten in the form 



m+1 1 m+1 "+ 



m+1 ^ m+1 

U_i 



M A 7 2 



, (B19) 



where a± = fx (s' ± A) are the eigenvalues of A. The L2 
solution of the homogeneous wave equation (|B6[) is given 



by (|B10p . Nevertheless, the integration constant a satis- 
fies a™ +1 u = g. It can be shown that the system (IB6|) 
with BCs (|B19[) is boundary stable and, according to [23|, 
the solution satisfies the following estimate 



/ni + 1 . , 1 1 b ~r a 

" °°m-l 

- / E h ro_i ^ 



m+1 



2=0 



< c 



3=0 



x=0 



(B20) 



for some strictly positive constant C > 0. Thus, by in- 
verting the Laplace-Fourier transformation we can esti- 
mate the Li norm of higher derivatives of the solution in 
terms of given data. 



Appendix C: Implementation of boundary 
conditions in spherical symmetry 

In this appendix we describe the numerical implemen- 
tation of the second order constraint preserving boundary 
conditions (|60M62p in spherical symmetry. 

We write the spherical line-element as 

ds 2 = X^lrr&T 2 + X-^TT-W, (CI) 

where dfl 2 = d8 2 + sin 2 9d(j) 2 . Similarly we evolve 
(K, A rr , At) for the extrinsic curvature. In spherical 
symmetry the algebraic constraints (|25p are 



D = log (> r 7 2 ) = 0, T=4^ 



A T 
2— = 0. 

7t 



(C2) 



Using the linearized equations of motion for the system, 
we rewrite the the boundary conditions as 



d r e--e, 

r 



2 r 



^^T^ V/3V 

- ~9 r 6 - ^d r K, 

d t k = - V2d r k- — k 

r 



4 

3r 2P 



' Or 



(C3) 

(C4) 
(C5) 



1 



6 - 

-A rr — drdrjT — ■^9,T r 

r 2 



- f (-2 + V2)9 r X + ^9 - 4(>r - 7t) 



drlT - d r \) 



2d r a--e~-—k~2T r 
3 3 



(C6) 
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where here, for brevity, we have linearized around flat 
space. The boundary condition for At can be obtained 
by using the algebraic constraints. Note the similarity 
with the Sommerfeld boundary condition. For the nu- 
merical implementation we populate ghostzones for each 
gridfunction /, by sixth order extrapolation 

/iV+i = 6/at+i-I — 15/tv+;-2 + 20/jV+i-3 — 15/AT+i-4 
+ QfN+i-5 ~ /iV+i-6) (C7) 



in order to approximate derivatives and compute Kreiss- 
Oliger artificial dissipation at the boundary. Here 
N denotes the boundary point. For the variables 
(0, T r , K, A rr , At) we simply replace the standard evo- 
lution equations with (|C3IIC6[) at the boundary. The re- 
maining variables are evolved according to their standard 
equation of motion at the boundary. 
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